Load preprocessed dataset (preprocessing code in 20_03_30_data_preprocessing.Rmd) and clustering (pipeline in 20_03_31_WGCNA.Rmd)
# Gandal dataset
load('./../Data/preprocessed_data.RData')
datExpr = datExpr %>% data.frame
DE_info = DE_info %>% data.frame
# GO Neuronal annotations: regex 'neuron' in GO functional annotations and label the genes that make a match as neuronal
GO_annotations = read.csv('./../Data/genes_GO_annotations.csv')
GO_neuronal = GO_annotations %>% filter(grepl('neuron', go_term)) %>%
mutate('ID'=as.character(ensembl_gene_id)) %>%
dplyr::select(-ensembl_gene_id) %>% distinct(ID) %>%
mutate('Neuronal'=1)
# SFARI Genes
SFARI_genes = read_csv('./../../../SFARI/Data/SFARI_genes_08-29-2019_w_ensembl_IDs.csv')
SFARI_genes = SFARI_genes[!duplicated(SFARI_genes$ID) & !is.na(SFARI_genes$ID),]
# Clusterings
clusterings = read_csv('./../Data/clusters.csv')
# Update DE_info with SFARI and Neuronal information
genes_info = DE_info %>% mutate('ID'=rownames(.)) %>% left_join(SFARI_genes, by='ID') %>%
mutate(`gene-score`=ifelse(is.na(`gene-score`), 'None', `gene-score`)) %>%
left_join(GO_neuronal, by='ID') %>% left_join(clusterings, by='ID') %>%
mutate(Neuronal=ifelse(is.na(Neuronal), 0, Neuronal)) %>%
mutate(gene.score=ifelse(`gene-score`=='None' & Neuronal==1, 'Neuronal', `gene-score`),
significant=padj<0.05 & !is.na(padj))
rm(DE_info, GO_annotations, clusterings)
print(paste0('Dynamic Tree leaves ', sum(genes_info$DynamicTree=='gray'), ' genes without cluster (',
round(mean(genes_info$DynamicTree=='gray')*100), '%)'))
## [1] "Dynamic Tree leaves 7841 genes without cluster (57%)"
print(paste0('Dynamic Hybrid leaves ', sum(genes_info$DynamicHybrid=='gray'), ' genes without cluster (',
round(mean(genes_info$DynamicHybrid=='gray')*100,2), '%)'))
## [1] "Dynamic Hybrid leaves 388 genes without cluster (2.83%)"
Dynamic Tree leaves more genes without a cluster, but in previous experiments it returned cleaner results, so I’m going to see which genes are lost to see how big the damage is
There seems to be a relation between DE and module membership, being DE a much more restrictive condition than being assigned to a cluster
pca = datExpr %>% prcomp
plot_data = data.frame('ID'=rownames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2]) %>%
left_join(genes_info, by='ID') %>% mutate('hasCluster'=DynamicTree!='gray',
'hasSFARIScore'=`gene-score`!='None') %>%
apply_labels(`gene-score`='SFARI Gene score', DynamicTree = 'Dynamic Tree Algorithm',
significant = 'Differentially Expressed', hasCluster = 'Belongs to a Module',
hasSFARIScore = 'Has a SFARI Score', syndromic = 'Has syndromic tag')
p1 = plot_data %>% ggplot(aes(PC1, PC2, color=hasCluster)) + geom_point(alpha=0.3) +
theme_minimal() + ggtitle('Genes are assigned to a cluster') + theme(legend.position='bottom')
p2 = plot_data %>% ggplot(aes(PC1, PC2, color=significant)) + geom_point(alpha=0.3) +
theme_minimal() + ggtitle('Genes found to be DE') + theme(legend.position='bottom')
grid.arrange(p1, p2, nrow=1)
rm(pca, p1, p2)
Almost all genes that don’t have a cluster are not differentially expressed
cat(paste0(round(100*sum(!plot_data$hasCluster & !plot_data$significant)/sum(!plot_data$hasCluster),1),
'% of the genes that don\'t have a cluster are not differentially expressed\n'))
## 98.1% of the genes that don't have a cluster are not differentially expressed
cro(plot_data$significant, list(plot_data$hasCluster, total()))
| Belongs to a Module | #Total | |||
|---|---|---|---|---|
| FALSE | TRUE | |||
| Differentially Expressed | ||||
| FALSE | 7690 | 5064 | 12754 | |
| TRUE | 151 | 793 | 944 | |
| #Total cases | 7841 | 5857 | 13698 | |
Most of the genes with a SFARI score are assigned to a cluster
cat(paste0(sum(plot_data$hasSFARIScore & !plot_data$hasCluster), ' of the SFARI genes (~',
round(100*sum(plot_data$hasSFARIScore & !plot_data$hasCluster)/sum(plot_data$hasSFARIScore)),
'%) are not assigned to any cluster\n'))
## 404 of the SFARI genes (~50%) are not assigned to any cluster
cro(plot_data$hasSFARIScore, list(plot_data$hasCluster, total()))
| Belongs to a Module | #Total | |||
|---|---|---|---|---|
| FALSE | TRUE | |||
| Has a SFARI Score | ||||
| FALSE | 7437 | 5455 | 12892 | |
| TRUE | 404 | 402 | 806 | |
| #Total cases | 7841 | 5857 | 13698 | |
The main difference between algorithms is that Dynamic Hybrid clusters outlier genes and Dynamic Tree leaves them out, so Dynamic Tree would give me a ‘cleaner’ group of genes to work with, without losing many SFARI genes, but Dynamic Hybrid has less and more balanced clusters
I think both options could be feasible, but I’m going to use the Dynamic Hybrid algorithm to keep more genes
clustering_selected = 'DynamicHybrid'
genes_info$Module = genes_info[,clustering_selected]
*The colour of the modules is the arbitrary one assigned during the WGCNA algorithm, where the gray cluster actually represents all the genes that were left without a cluster (so it’s not actually a cluster).
cat(paste0('The Dynamic Hybrid algorithm created ', length(unique(genes_info$Module))-1, ' modules and leaves ',
sum(genes_info$Module=='gray'), ' genes without a module.\n'))
## The Dynamic Hybrid algorithm created 102 modules and leaves 388 genes without a module.
sort(table(genes_info$Module), decreasing = TRUE)
##
## #45B500 #D277FF #00C0BE #CC9600 #29A3FF #00BF7F #00C0B8 gray #00C087 #FF6D8F
## 1403 833 774 640 485 435 413 388 333 332
## #89AC00 #FE61CD #FF61C7 #619CFF #00BAE0 #E08A00 #9091FF #80AD00 #E7851E #FF689E
## 305 302 266 199 194 183 177 172 155 154
## #FA62D9 #00BB46 #CB7AFF #FF65AD #30B600 #BA82FF #C79800 #F166E8 #E96AF1 #00BF77
## 152 148 146 143 139 139 139 137 131 125
## #98A800 #77AF00 #00B5EE #00BECA #FC717F #7398FF #F863DE #00BD5C #00B927 #62B200
## 123 120 115 115 113 109 108 104 96 96
## #EE8042 #D49200 #E48800 #D09400 #B79F00 #00B709 #FC61D3 #FF62BA #FF6A97 #00A6FF
## 91 90 90 87 84 81 74 73 73 72
## #A5A500 #00C1A4 #A789FF #DF70F9 #F564E3 #FA7476 #C29B00 #55B300 #D973FC #EB8333
## 71 68 68 68 65 65 64 63 57 57
## #00B0F6 #9FA700 #ED68ED #00BD65 #00C1AB #00A9FF #BD9D00 #F07E4E #00AEF9 #8295FF
## 55 54 53 52 51 50 50 50 49 49
## #00BCD6 #ACA300 #FD6F87 #00B8E5 #FF63B4 #00C19D #9C8DFF #F37B59 #00C08E #FF67A6
## 48 47 45 43 43 42 41 41 40 40
## #00B3F2 #00ABFC #00BE6E #00C0B2 #FF61C1 #4B9FFF #00BC51 #00BFC4 #90AA00 #6DB100
## 39 38 36 36 36 34 33 33 33 32
## #00B7E9 #DD8D00 #F67963 #00BBDB #B186FF #D98F00 #00C096 #E46DF5 #B2A100 #C37EFF
## 30 30 28 27 25 25 24 24 22 22
## #F8766D #00BDD0 #00BA38
## 18 16 12
plot_data = table(genes_info$Module) %>% data.frame %>% arrange(desc(Freq))
ggplotly(plot_data %>% ggplot(aes(x=reorder(Var1, -Freq), y=Freq)) + geom_bar(stat='identity', fill=plot_data$Var1) +
ggtitle('Module size') + ylab('Number of genes') + xlab('Module') + theme_minimal() +
theme(axis.text.x = element_text(angle = 90)))
In the WGCNA documentation they use Pearson correlation to calculate correlations, I think all of their variables were continuous. Since I have categorical variables I’m going to use the hetcor function, that calculates Pearson, polyserial or polychoric correlations depending on the type of variables involved.
I’m not sure how the corPvalueStudent function calculates the p-values and I cannot find any documentation…
Compared correlations using Pearson correlation and with hetcor and they are very similar, but a bit more extreme with hetcor. The same thing happens with the p-values.
datTraits = datMeta %>% dplyr::select(Diagnosis, Gender, Age, brain_lobe, PMI, SiteHM) %>%
dplyr::rename('ExtractionBatch' = SiteHM)
# Recalculate MEs with color labels
ME_object = datExpr %>% t %>% moduleEigengenes(colors = genes_info$Module)
MEs = orderMEs(ME_object$eigengenes)
# Calculate correlation between eigengenes and the traits and their p-values
moduleTraitCor = MEs %>% apply(2, function(x) hetcor(x, datTraits)$correlations[1,-1]) %>% t
rownames(moduleTraitCor) = colnames(MEs)
colnames(moduleTraitCor) = colnames(datTraits)
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nrow(datExpr))
# Create text matrix for the Heatmap
textMatrix = paste0(signif(moduleTraitCor, 2), ' (', signif(moduleTraitPvalue, 1), ')')
dim(textMatrix) = dim(moduleTraitCor)
# In case there are any NAs
if(sum(!complete.cases(moduleTraitCor))>0){
print(paste0(sum(is.na(moduleTraitCor)),' correlation(s) could not be calculated'))
}
rm(ME_object)
Note: In Gandal’s dataset there were some problems calculating the Module-Diagnosis correlation for Modules with a correlation too strong with Diagnosis, this doesn’t seem to be a problem here, probably because the dataset is a lot noisier and the Diagnosis groups don’t separate as clearly as with the Gandal dataset (as we saw in the Sample’s PCA plot in the exploratory analysis)
# Calculate the correlation tha failed with hetcor()
missing_modules = rownames(moduleTraitCor)[is.na(moduleTraitCor[,1])]
for(m in missing_modules){
cat(paste0('Correcting Module-Diagnosis correlation for Module ', m))
moduleTraitCor[m,'Diagnosis'] = polyserial(MEs[,m], datTraits$Diagnosis)
}
rm(missing_modules)
Modules have very strong correlations with Diagnosis with really small p-values and not much relation with anything else. Perhaps a little with Batch, probably because we didn’t do the batch correction because the batch was correlated to Diagnosis and we didn’t want to remove any biological signal
They gray ‘module’ has one of the lowest correlations with Diagnosis
# Sort moduleTraitCor by Diagnosis
moduleTraitCor = moduleTraitCor[order(moduleTraitCor[,1], decreasing=TRUE),]
moduleTraitPvalue = moduleTraitPvalue[order(moduleTraitCor[,1], decreasing=TRUE),]
# Create text matrix for the Heatmap
textMatrix = paste0(signif(moduleTraitCor, 2), ' (', signif(moduleTraitPvalue, 1), ')')
dim(textMatrix) = dim(moduleTraitCor)
labeledHeatmap(Matrix = moduleTraitCor, xLabels = names(datTraits), yLabels = gsub('ME','',rownames(moduleTraitCor)),
yColorWidth=0, colors = brewer.pal(11,'PiYG'), bg.lab.y = gsub('ME','',rownames(moduleTraitCor)),
textMatrix = textMatrix, setStdMargins = FALSE, cex.text = 0.8, cex.lab.y = 0.75, zlim = c(-1,1),
main = paste('Module-Trait relationships'))
diagnosis_cor = data.frame('Module' = gsub('ME','',rownames(moduleTraitCor)),
'MTcor' = moduleTraitCor[,'Diagnosis'],
'MTpval' = moduleTraitPvalue[,'Diagnosis'])
genes_info = genes_info %>% left_join(diagnosis_cor, by='Module')
rm(moduleTraitCor, moduleTraitPvalue, textMatrix, diagnosis_cor)
Modules with a high Module-Diagnosis (absolute) correlation have a higher content of differentially expressed genes
plot_data = genes_info %>% group_by(Module, MTcor) %>% summarise(p = 100*mean(significant))
plot_data %>% ggplot(aes(MTcor, p)) + geom_hline(yintercept=mean(plot_data$p), color='gray', linetype='dotted') +
geom_point(color=plot_data$Module, aes(id=Module)) + theme_minimal() +
xlab('Modules ordered by Module-Diagnosis correlation') + ylab('Percentage of differentially expressed genes')
Gene significance: is the value between the correlation between the gene and the trait we are interested in. A positive gene significance means the gene is overexpressed and a negative value means its underexpressed. (The term ‘significance’ is not very acurate because it’s not actually measuring statistical significance, it’s just a correlation, but that’s how they call it in WGCNA…)
Module Membership is the correlation of the module’s eigengene and the expression profile of a gene. The higher the Module Membership, the more similar the gene is to the genes that constitute the module. (I won’t use this measure yet)
# It's more efficient to iterate the correlations one by one, otherwise it calculates correlations between the eigengenes and also between the genes, which we don't need
# Check if MM information already exists and if not, calculate it
if(file.exists(paste0('./../Data/dataset_', clustering_selected, '.csv'))){
dataset = read.csv(paste0('./../Data/dataset_', clustering_selected, '.csv'))
dataset$Module = dataset[,clustering_selected]
} else {
############# 1. Calculate Gene Significance
GS_info = data.frame('ID' = rownames(datExpr),
'GS' = datExpr %>% apply(1, function(x) hetcor(x, datMeta$Diagnosis)$correlations[1,2])) %>%
mutate('GSpval' = corPvalueStudent(GS, ncol(datExpr)))
############# 2. Calculate Module Membership
#setup parallel backend to use many processors
cores = detectCores()
cl = makeCluster(cores-1)
registerDoParallel(cl)
# Create matrix with MM by gene
MM = foreach(i=1:nrow(datExpr), .combine=rbind) %dopar% {
library(polycor)
tempMatrix = apply(MEs, 2, function(x) hetcor(as.numeric(datExpr[i,]), x)$correlations[1,2])
tempMatrix
}
# Stop clusters
stopCluster(cl)
rownames(MM) = rownames(datExpr)
colnames(MM) = paste0('MM',gsub('ME','',colnames(MEs)))
# Calculate p-values
MMpval = MM %>% corPvalueStudent(ncol(datExpr)) %>% as.data.frame
colnames(MMpval) = paste0('MMpval', gsub('ME','',colnames(MEs)))
MM = MM %>% as.data.frame %>% mutate(ID = rownames(.))
MMpval = MMpval %>% as.data.frame %>% mutate(ID = rownames(.))
# Join and save results
dataset = genes_info %>% dplyr::select(ID, `gene-score`, clustering_selected, MTcor, MTpval) %>%
left_join(GS_info, by='ID') %>%
left_join(MM, by='ID') %>%
left_join(MMpval, by='ID')
write.csv(dataset, file = paste0('./../Data/dataset_', clustering_selected, '.csv'), row.names = FALSE)
rm(cores, cl)
}
GS_missing = dataset$ID[is.na(dataset$GS)] %>% as.character
if(length(GS_missing)>0){
print(paste0(length(GS_missing),' correlations between genes and Diagnosis could not be calculated, ',
'calculating them with the polyserial function'))
for(g in GS_missing){
dataset$GS[dataset$ID == g] = polyserial(as.numeric(datExpr[g,]), datMeta$Diagnosis)
}
}
rm(GS_missing)
Gene significance and Log Fold Chance are two different ways to measure the same thing, so there should be a concordance between them
plot_data = dataset %>% dplyr::select(ID, MTcor, GS) %>% left_join(genes_info %>% dplyr::select(ID, gene.score), by='ID') %>%
left_join(genes_info %>% dplyr::select(ID, baseMean, log2FoldChange, significant, Module), by='ID') %>%
left_join(data.frame(MTcor=unique(dataset$MTcor)) %>% arrange(by=MTcor) %>%
mutate(order=1:length(unique(dataset$MTcor))), by='MTcor')
ggplotly(plot_data %>% ggplot(aes(GS, log2FoldChange)) + geom_point(color=plot_data$Module, alpha=0.5, aes(ID=Module)) +
geom_smooth(color='gray') + theme_minimal() + xlab('Gene Significance') +
ggtitle(paste0('Correlation = ', round(cor(plot_data$log2FoldChange, plot_data$GS)[1], 4))))
In general, modules with the highest Module-Diagnosis correlation should have genes with high Gene Significance
Note: For the Module-Diagnosis plots, if you do boxplots, you lose the exact module-diagnosis correlation and you only keep the order, so I decided to compensate this downside with a second plot, where each point is plotted individually using their module’s Module-Diagnosis correlation as the x axis. I think the boxplot plot is easier to understand but the second plot contains more information, so I don’t know which one is better.
plot_data = plot_data %>% arrange(order)
ggplotly(plot_data %>% ggplot(aes(order, GS, group=order)) + geom_hline(yintercept=0, color='gray', linetype='dotted') +
geom_boxplot(fill=unique(plot_data$Module)) + theme_minimal() +
xlab('Modules ordered by Module-Diagnosis correlation') + ylab('Gene Significance'))
plot_data %>% ggplot(aes(MTcor, GS)) + geom_hline(yintercept=0, color='gray', linetype='dotted') +
geom_point(color=plot_data$Module, alpha=0.1, aes(id=ID)) + geom_smooth(color='gray', alpha=0.3) +
theme_minimal() + xlab('Module-Diagnosis correlation') + ylab('Gene Significance') +
ggtitle(paste0('R^2=',round(cor(plot_data$MTcor, plot_data$GS)^2,4)))
The same should happen with the Log Fold Change
ggplotly(plot_data %>% ggplot(aes(order, log2FoldChange, group=order)) + geom_hline(yintercept=0, color='gray', linetype='dotted') +
geom_boxplot(fill=unique(plot_data$Module)) +
theme_minimal() + xlab('Modules ordered by Module-Diagnosis correlation') + ylab('log2FoldChange'))
ggplotly(plot_data %>% ggplot(aes(MTcor, log2FoldChange)) + geom_hline(yintercept=0, color='gray', linetype='dotted') +
geom_point(color=plot_data$Module, alpha=0.1, aes(id=ID)) + geom_smooth(color='gray', alpha=0.3) +
theme_minimal() + xlab('Module-Diagnosis correlation') + ylab('log2FoldChange') +
ggtitle(paste0('R^2=',round(cor(plot_data$MTcor, plot_data$log2FoldChange)^2,4))))
We can see a small relation between module-Diagnosis and mean expression that could be explained by what we had observed in 20_03_31_exploratory_analysis.html where overexpressed genes tended to have lower levels of expression than the overexpressed genes, but the pattern is not very strong
ggplotly(plot_data %>% ggplot(aes(order, log2(baseMean+1), group=order)) +
geom_hline(yintercept=mean(log2(plot_data$baseMean+1)), color='gray', linetype='dotted') +
geom_boxplot(fill=unique(plot_data$Module)) + theme_minimal() +
xlab('Modules ordered by Module-Diagnosis correlation') + ylab('log2(Mean Expression)'))
plot_data %>% ggplot(aes(MTcor, log2(baseMean+1))) + geom_point(alpha=0.2, color=plot_data$Module, aes(id=ID)) +
geom_hline(yintercept=mean(log2(plot_data$baseMean+1)), color='gray', linetype='dotted') +
geom_smooth(color='gray', alpha=0.3) + theme_minimal() + xlab('Module-Diagnosis correlation') +
ggtitle(paste0('R^2=',round(cor(plot_data$MTcor, log2(plot_data$baseMean+1))^2,4)))
All of the variables seem to agree with each other, Modules with a high correlation with Diagnosis tend to have genes with high values of Log Fold Change as well as high values of Gene Significance, and the gray module, which groups all the genes that weren’t assigned to any cluster tends to have a very poor performance in all of the metrics.
Since SFARI scores genes depending on the strength of the evidence linking it to the development of autism, in theory, there should be some concordance between the metrics we have been studying above and these scores…
Genes with SFARI score 6 have the highest Gene Significance median of all groups (!)
There doesn’t seem to be a big difference between Gene significance distributions and SFARI scores, or even Neuronal tag and the rest of the genes
There isn’t a clear pattern within the SFARI genes as the one we could see in Gandal’s dataset
ggplotly(plot_data %>% ggplot(aes(gene.score, abs(GS), fill=gene.score)) + geom_boxplot() +
scale_fill_manual(values=SFARI_colour_hue(r=c(1:6,8,7))) + theme_minimal() +
ylab('abs(Gene Significance)') + xlab('SFARI Scores') + theme(legend.position='none'))
All groups have very similar distributions with relatively high Module-Diagnosis correlations
There could be a negative pattern between SFARI score and Module-diagnosis correlation, but it’s not as clear as the one we could see in Gandal’s dataset
ggplotly(plot_data %>% ggplot(aes(gene.score, abs(MTcor), fill=gene.score)) + geom_boxplot() +
scale_fill_manual(values=SFARI_colour_hue(r=c(1:6,8,7))) + theme_minimal() +
ylab('abs(Module-Diagnosis Correlation)') + xlab('SFARI Scores') + theme(legend.position='none'))
SFARI genes are not very consistent with the other measurements, but they don’t seem to contradict them in such as strong way as we had seen in Gandal’s dataset. Here it looks like they may be relatively independent from each other.
sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
##
## locale:
## [1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB.UTF-8 LC_COLLATE=en_GB.UTF-8
## [5] LC_MONETARY=en_GB.UTF-8 LC_MESSAGES=en_GB.UTF-8
## [7] LC_PAPER=en_GB.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] doParallel_1.0.15 iterators_1.0.12 foreach_1.5.0
## [4] polycor_0.7-10 expss_0.10.2 WGCNA_1.69
## [7] fastcluster_1.1.25 dynamicTreeCut_1.63-1 GGally_1.5.0
## [10] gridExtra_2.3 viridis_0.5.1 viridisLite_0.3.0
## [13] RColorBrewer_1.1-2 dendextend_1.13.4 plotly_4.9.2
## [16] glue_1.3.2 reshape2_1.4.3 forcats_0.5.0
## [19] stringr_1.4.0 dplyr_0.8.5 purrr_0.3.3
## [22] readr_1.3.1 tidyr_1.0.2 tibble_3.0.0
## [25] ggplot2_3.3.0 tidyverse_1.3.0
##
## loaded via a namespace (and not attached):
## [1] colorspace_1.4-1 ellipsis_0.3.0
## [3] htmlTable_1.13.3 XVector_0.24.0
## [5] GenomicRanges_1.36.1 base64enc_0.1-3
## [7] fs_1.3.2 rstudioapi_0.11
## [9] farver_2.0.3 bit64_0.9-7
## [11] mvtnorm_1.1-0 AnnotationDbi_1.46.1
## [13] fansi_0.4.1 lubridate_1.7.4
## [15] xml2_1.2.5 codetools_0.2-16
## [17] splines_3.6.3 impute_1.58.0
## [19] geneplotter_1.62.0 knitr_1.28
## [21] Formula_1.2-3 jsonlite_1.6.1
## [23] annotate_1.62.0 broom_0.5.5
## [25] cluster_2.1.0 GO.db_3.8.2
## [27] dbplyr_1.4.2 png_0.1-7
## [29] compiler_3.6.3 httr_1.4.1
## [31] backports_1.1.5 assertthat_0.2.1
## [33] Matrix_1.2-18 lazyeval_0.2.2
## [35] cli_2.0.2 acepack_1.4.1
## [37] htmltools_0.4.0 tools_3.6.3
## [39] GenomeInfoDbData_1.2.1 gtable_0.3.0
## [41] Rcpp_1.0.4 Biobase_2.44.0
## [43] cellranger_1.1.0 vctrs_0.2.4
## [45] preprocessCore_1.46.0 nlme_3.1-144
## [47] crosstalk_1.1.0.1 xfun_0.12
## [49] rvest_0.3.5 lifecycle_0.2.0
## [51] XML_3.99-0.3 zlibbioc_1.30.0
## [53] scales_1.1.0 hms_0.5.3
## [55] SummarizedExperiment_1.14.1 yaml_2.2.1
## [57] memoise_1.1.0 rpart_4.1-15
## [59] reshape_0.8.8 latticeExtra_0.6-29
## [61] stringi_1.4.6 RSQLite_2.2.0
## [63] genefilter_1.66.0 S4Vectors_0.22.1
## [65] checkmate_2.0.0 BiocGenerics_0.30.0
## [67] BiocParallel_1.18.1 GenomeInfoDb_1.20.0
## [69] bitops_1.0-6 rlang_0.4.5
## [71] pkgconfig_2.0.3 matrixStats_0.56.0
## [73] evaluate_0.14 lattice_0.20-40
## [75] labeling_0.3 htmlwidgets_1.5.1
## [77] bit_1.1-15.2 tidyselect_1.0.0
## [79] plyr_1.8.6 magrittr_1.5
## [81] DESeq2_1.24.0 R6_2.4.1
## [83] IRanges_2.18.3 generics_0.0.2
## [85] Hmisc_4.4-0 DelayedArray_0.10.0
## [87] DBI_1.1.0 mgcv_1.8-31
## [89] pillar_1.4.3 haven_2.2.0
## [91] foreign_0.8-75 withr_2.1.2
## [93] RCurl_1.98-1.1 survival_3.1-11
## [95] nnet_7.3-13 modelr_0.1.6
## [97] crayon_1.3.4 rmarkdown_2.1
## [99] jpeg_0.1-8.1 locfit_1.5-9.4
## [101] grid_3.6.3 readxl_1.3.1
## [103] data.table_1.12.8 blob_1.2.1
## [105] reprex_0.3.0 digest_0.6.25
## [107] xtable_1.8-4 stats4_3.6.3
## [109] munsell_0.5.0